function nhfile = mtc_Resample(hfile, ifunc, ts)
% MTC::Resample  - temporally resample an MTC
%
% FORMAT:       newmtc = mtc.Resample(ifunc, ts)
%
% Input fields:
%
%       ifunc       interpolation function,
%                   'linear', 'cubic', 'nearest', 'lanczos3'
%       ts          temporal sampling points (in MS!)
%
% Output fields:
%
%       newmtc      MTC with resampled timepoints

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'mtc') || ...
   ~ischar(ifunc) || ...
   ~any(strcmpi(ifunc(:)', {'cubic', 'lanczos3', 'linear', 'nearest'})) || ...
   ~isa(ts, 'double') || ...
    isempty(ts) || ...
    length(ts) ~= numel(ts) || ...
    (numel(ts) > 2 && any(abs(diff(ts, 2) > 1))) || ...
    any(isinf(ts(:)) | isnan(ts(:)))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
ifunc = lower(ifunc(:)');
bc = bvqxfile_getcont(hfile.L);
nv = size(bc.MTCData, 1);
np = size(bc.MTCData, 2);
if any(ts < 0 | ts > (bc.TR * (nv - 0.5)))
    warning( ...
        'BVQXfile:BadArgument', ...
        'Sampling points out of bounds.' ...
    );
end

% get newTR
if numel(ts) > 1
    nTR = round(mean(diff(ts(:))));
else
    nTR = bc.TR;
end
TRfac = nTR / bc.TR;
ts = 1 + ts(:) / bc.TR;
sts = min(np, 4096);
tss = [Inf, Inf; ts(1), 1; mean(diff(ts)), 1; ts(end), np];

% try progress bar
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', 'Resampling MTC...');
    BVQXprogress(pbar, 0, 'Preparing structures...', 'visible', 0, 10);
    pvx = 9 / np;
catch
    pbar = [];
end

% create data arrays (old in single AND new in uint16 must fit into mem!)
nMTC = single(0);
try
    oMTC = bc.MTCData(:, :);
    nMTC(numel(ts), np) = nMTC(1);
catch
    error( ...
        'BVQXfile:OutOfMemory', ...
        'Out of memory (old and new VTC must fit in memory).' ...
    );
end

% create new VTC object
if ~isempty(pbar)
    BVQXprogress(pbar, 0.5, 'Creating new object...');
end
nhfile = BVQXfile('new:mtc');
ncfile = bvqxfile_getcont(nhfile.L);
ncfile.NrOfVertices = size(bc.MTCData, 2);
ncfile.NrOfTimePoints = numel(ts);
ncfile.SourceVTCFile = bc.SourceVTCFile;
ncfile.LinkedPRTFile = bc.LinkedPRTFile;
ncfile.HemodynamicDelay = bc.HemodynamicDelay;
ncfile.TR = nTR;
ncfile.HRFDelta = bc.HRFDelta;
ncfile.HRFTau = bc.HRFTau;
ncfile.ProtocolSegmentSize = fix(bc.ProtocolSegmentSize * TRfac);
ncfile.ProtocolSegmentOffset = fix(bc.ProtocolSegmentOffset * TRfac);
ncfile.MTCData = single([]);

% loop over points with stepsize
if ~isempty(pbar)
    BVQXprogress(pbar, 1, sprintf('Resampling with %s interpolation...', ifunc));
end
for sc = 1:sts:np
    tss(2, 2) = sc;
    tss(4, 2) = min(np, sc + sts - 1);
    nMTC(:, tss(2, 2):tss(4, 2)) = ...
        single(flexinterpn_method(oMTC, tss, 0, ifunc));
    if ~isempty(pbar)
        BVQXprogress(pbar, 1 + sc * pvx);
    end
end

% set final field
ncfile.MTCData = nMTC;
clear oMTC;

% but content into array
bvqxfile_setcont(nhfile.L, ncfile);

% remove bar
if ~isempty(pbar)
    closebar(pbar);
end
